# ----------------------------------------------------------------------
# How Does Shaming Human Rights Violators Abroad Shape Attitudes at Home?
# Study II: Script to produce figures for manuscript and SM
# ----------------------------------------------------------------------

# ----------------------------------------------------------------------
# load all relevant themes and packages
# ----------------------------------------------------------------------
rm(list=ls())
library("tidyverse")
library("ggplot2")
library("stargazer")
library("xtable")
library("texreg")
library("reshape2")
library("RItools")
library("estimatr")
library("fastDummies")
library("ggpubr")
library("estimatr")
library("mediation")
library("car")
library("table1")
library("psy")
library("dplyr")
library("ggstatsplot")
library("ggplot2")
library("stringr")
library("tidyverse")
library("haven")
library("lfe")
library("naniar")
library("fixest")
library("Hmisc")
library("quanteda")
library("readtext")
library("htmltools")
library("devtools")
library("quanteda.textplots")

# ----------------------------------------------------------------------
# Read data
# ----------------------------------------------------------------------

data <- read_csv("Data_S2.csv")
data <- data[3:nrow(data),]

# ----------------------------------------------------------------------
# Re-code covariates and outcomes
# ----------------------------------------------------------------------

#OUTCONES summary and rescaling

#`gov_approve` -- rescale: high # indicates support for government
data$gov_approve = recode(data$gov_approve1, "1=7; 2=6; 3=5; 4=4; 5=3; 6=2; 7=1")

#`torture` -- high # indicate opposition to torture (no need to rescale)

# `legal` -- index: create index of legal obligation measure (Bayram 2017) 

data <- data %>% 
  mutate(.,
         legal = (as.numeric(legal_obligation_1) + 
                    as.numeric(legal_obligation_2) + 
                    as.numeric(legal_obligation_3) + 
                    as.numeric(legal_obligation_4) + 
                    as.numeric(legal_obligation_5))/5)
data$legal<-6-data$legal

# `HR values` -- high # indicates more respect for human rights (no change required)

# morality -- rescale: high # indicates U.S. more moral 
data$morality = recode(data$morality, "1=7; 2=6; 3=5; 4=4; 5=3; 6=2; 7=1")

# power -- rescale: high # indicates U.S. more powerful 
data$power = recode(data$power, "1=7; 2=6; 3=5; 4=4; 5=3; 6=2; 7=1")

# ----------------------------------------------------------------------
# recode variables
# ----------------------------------------------------------------------

data <- data %>%
  mutate(.,
         treat_shame = case_when(
           condition == "no_shame_ally" ~ 0,
           condition == "no_shame_no_ally" ~ 0,
           condition == "shame_ally" ~ 1,
           condition == "shame_no_ally" ~ 1),
         treat_target = case_when(
           condition == "no_shame_ally" ~ 1,
           condition == "no_shame_no_ally" ~ 0,
           condition == "shame_ally" ~ 1,
           condition == "shame_no_ally" ~ 0),
      Democrat = case_when(
           partisanship == 1 ~ 0,
           partisanship == 2 ~ 1,
           partisanship == 3 ~ 0),
         Republican = case_when(
           partisanship == 1 ~ 1,
           partisanship == 2 ~ 0,
           partisanship == 3 ~ 0),
         Independents = case_when(
           partisanship == 1 ~ 0,
           partisanship == 2 ~ 0,
           partisanship == 3 ~ 1
         ),
         Female = ifelse(gender == 2, 1, 0),
         White = ifelse(race == 1, 1, 0),
         Black = ifelse(race == 2,1,0),
         Hispanic = ifelse(race==3,1,0),
         Asian = ifelse(race==4,1,0),
         Native_american = ifelse(race==5,1,0),
         Middel_eastern = ifelse(race==6,1,0),
         Mixed = ifelse(race==7,1,0))

# ----------------------------------------------------------------------
# Produce Figure 6 (manuscript)
# ----------------------------------------------------------------------

#create sub datas for ally and non-ally conditions
target_ally <- data[ which(data$treat_target==1), ]
target_nonally <- data[ which(data$treat_target==0), ]


#H1: does shaming affect government approval?
H1ally<- lm_robust(gov_approve~treat_shame , data = target_ally)%>% tidy %>% 
  mutate(.,
         model = "Government \n Approval \n (H1)",
         Target = "Ally") 

H1nonally<- lm_robust(gov_approve~treat_shame , data = target_nonally)%>% tidy %>% 
  mutate(.,
         model = "Government \n Approval \n (H1)",
         Target = "Non-Ally") 

#H2: does shaming affect human rights perception?
H2ally<- lm_robust(`HR values`~treat_shame, data=target_ally)%>% tidy %>% 
  mutate(.,
         model = "Human Rights \n Perception \n (H2)",
         Target = "Ally")

H2nonally<- lm_robust(`HR values`~treat_shame, data=target_nonally)%>% tidy %>% 
  mutate(.,
         model = "Human Rights \n Perception \n (H2)",
         Target = "Non-Ally")

#H3a: does shaming affect opposition of torture?
H3_ally<- lm_robust(torture~treat_shame, data = target_ally) %>% tidy %>% 
  mutate(.,
         model = "Oppose \n Torture \n (H3)",
         Target = "Ally")

H3_nonally<- lm_robust(torture~treat_shame, data = target_nonally) %>% tidy %>% 
  mutate(.,
         model = "Oppose \n Torture \n (H3)",
         Target = "Non-Ally")

#Create figure for main effects

main_fx <- rbind(H1ally, H1nonally, H2ally, H2nonally, 
                 H3_ally, H3_nonally) %>% 
  filter(.,
         term == "treat_shame")


order_x <- c("Government \n Approval \n (H1)", "Human Rights \n Perception \n (H2)",
             "Oppose \n Torture \n (H3)")

# Generate coefficient plot 
ggplot(main_fx, aes(x = model, y = estimate, color = Target, shape = Target)) +
  geom_hline(yintercept = 0, color = "gray50", linetype = 2, size = 0.2) +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high),
                  position = position_dodge(width = 0.4)) +
  #coord_flip()+
  labs(x = "",
       y = "Effect Size") +
  scale_x_discrete(limits= order_x) + 
  scale_color_manual(values = c("firebrick2","dodgerblue2")) +
  scale_fill_manual(values = c("firebrick2","dodgerblue2")) +
  theme(text = element_text(size = 10, family = "Times"),
        legend.key=element_blank(),
        panel.grid.major = element_blank(), 
        axis.text.x = element_text(size = 10),
        plot.caption = element_text(size = 10, family = "Times",hjust = -.02),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(), 
        axis.line = element_line(colour = "black"))

ggsave("../Main_Figs/fig6.pdf", width = 6, height = 4)


#substantive discussion of Fig 6:
#effect in terms of SDs:
1.4/sd(data$gov_approve,na.rm = T)
0.99/sd(data$`HR values`,na.rm = T)
0.2/sd(data$torture,na.rm=T)

# ----------------------------------------------------------------------
# Supplementary Materials
# ----------------------------------------------------------------------

#All files are saved into the Appendix/S2 folders

# ----------------------------------------------------------------------
# Descriptive statistics
# ----------------------------------------------------------------------

data1 <- data %>% 
  filter(.,
         !is.na(partisanship))

data1$age2<-as.numeric(data1$age2)
data1$gov_approve<-as.numeric(data1$gov_approve)
data1$`HR values`<-as.numeric(data1$`HR values`)
data1$torture<-as.numeric(data1$torture)
data1$legal<-as.numeric(data1$legal)
data1$power<-as.numeric(data1$power)
data1$morality<-as.numeric(data1$morality)

disc_table <- dplyr::select(data1, c(gov_approve,`HR values`, torture, legal, morality, power,
                                     Democrat, Independents, Republican, Female, White, Black,
                                     Hispanic, Asian, Native_american, Middel_eastern, Mixed, age2))

stargazer(as.data.frame(disc_table),
          covariate.labels = c("Government Approval","Human Rights Respect","Oppose Torture", "Legal Obligation", "Morality", "Power",
                               "Democrat", "Independent", "Republican", "Female", 
                               "White", "Black/African American", "Hispanic", "Asian/Asian American", "Native American",
                               "Middle Eastern", "Mixed Race","Age"),
          
          title = "Descriptive Statistics - Survey Experiment II",
          label = "tab:dsc_expII",
          style = "qje",
          notes.append = T,
          notes.align = "l",
          out="../Appendix/Tables/S2/B6.tex")


# ----------------------------------------------------------------------
# Tables for main models
# ----------------------------------------------------------------------

#Create table of main models and taeget and demographic controls
#Add (pre-registered) legal obligation outcome

H1<- lm(gov_approve~treat_shame, data = data)
H2<- lm(`HR values`~treat_shame, data=data)
H3a<- lm(torture~treat_shame, data = data)
H3b<- lm(legal~treat_shame, data = data)

##add target identity control

H1target<- lm(gov_approve~treat_shame+treat_target, data = data)
summary(H1target)

H2target<- lm(`HR values`~treat_shame+treat_target, data=data)
summary(H2target)

H3atarget<- lm(torture~treat_shame+treat_target, data = data)
summary(H3atarget)

H3btarget<- lm(legal~treat_shame+treat_target, data = data)
summary(H3btarget)

##add demographic controls

#Govermnemt approval
H1_control<- lm(gov_approve~treat_shame+treat_target + as.factor(gender) + as.factor(age2) 
                + as.factor(ethnicity) + as.factor(education) + as.factor(hispanic) +
                  as.factor(political_party) + as.factor(region), data = data) 
summary(H1_control)

#HR perception
H2_control <- lm(`HR values`~treat_shame+treat_target + as.factor(gender) + as.factor(age2) 
                 + as.factor(ethnicity) + as.factor(education) + as.factor(hispanic) +
                   as.factor(political_party) + as.factor(region), data = data) 
summary(H2_control)

#Torture
H3_conta<- lm(torture~treat_shame+treat_target + as.factor(gender) + as.factor(age2) 
              + as.factor(ethnicity) + as.factor(education) + as.factor(hispanic) +
                as.factor(political_party) + as.factor(region), data = data) 
summary(H3_conta)

#Legal
H3_contb<- lm(legal~treat_shame+treat_target + as.factor(gender) + as.factor(age2) 
              + as.factor(ethnicity) + as.factor(education) + as.factor(hispanic) +
                as.factor(political_party) + as.factor(region), data = data) 
summary(H3_contb)

##regress target*shaming + controls
#Govermnemt approval
H1_int<- lm(gov_approve~treat_shame*treat_target + as.factor(gender) + as.factor(age2) 
            + as.factor(ethnicity) + as.factor(education) + as.factor(hispanic) +
              as.factor(political_party) + as.factor(region), data = data) 
summary(H1_int)

#HR perception
H2_int <- lm(`HR values`~treat_shame*treat_target + as.factor(gender) + as.factor(age2) 
             + as.factor(ethnicity) + as.factor(education) + as.factor(hispanic) +
               as.factor(political_party) + as.factor(region), data = data) 
summary(H2_int)

#Torture
H3a_int<- lm(torture~treat_shame*treat_target + as.factor(gender) + as.factor(age2) 
             + as.factor(ethnicity) + as.factor(education) + as.factor(hispanic) +
               as.factor(political_party) + as.factor(region), data = data) 
summary(H3a_int)

#Legal
H3b_int<- lm(legal~treat_shame*treat_target + as.factor(gender) + as.factor(age2) 
             + as.factor(ethnicity) + as.factor(education) + as.factor(hispanic) +
               as.factor(political_party) + as.factor(region), data = data) 
summary(H3b_int)


##Create tables

#Create table for H1
tab1<-stargazer(H1, H1target, H1_control, H1_int,
                out= "../Appendix/Tables/S2/B7.tex",  style="qje",
                title="Treatment effects on government approval (H1)",
                keep = c("treat_shame", "treat_target", "treat_shame*treat_target"),
                dep.var.labels = "Government Approval",
                covariate.labels = c("Shaming treatment", "Ally treatment", "Shaming*Ally"),
                omit.stat = c("rsq", "f", "ser", "adj.rsq"),
                star.cutoffs = NULL,
                label = "tab:H1models_II",
               # star.char = c("", "", ""),
                add.lines = list(
                  c("Demographic Controls",  "No", "No", "Yes", "Yes")),
                ci=FALSE,
                omit.table.layout = NULL,
                notes.append = FALSE, 
                notes.align = "l",
                notes = c(""))

note.latex <- "\\multicolumn{5}{c} {\\parbox[t]{12cm}{ \\textit{Notes:} Demographic controls include: sex, age, race,
ethnicity, education, ideology, and region.
* $p < .1$, ** $p < .05$, *** $p < .01$}} \\\\"
tab1[grepl("Note",tab1)] <- note.latex
cat(tab1, sep = "\n")
write(tab1, '../Appendix/Tables/S2/B7.tex')

#Create table for H2
tab1<-stargazer(H2, H2target, H2_control, H2_int,
                out= "../Appendix/Tables/S2/B8.tex",  style="qje",
                title="Treatment effects on human rights perceptions (H2)",
                keep = c("treat_shame", "treat_target", "treat_shame*treat_target"),
                dep.var.labels = "Human Rights Perceptions",
                covariate.labels = c("Shaming treatment", "Ally treatment", "Shaming*Ally"),
                omit.stat = c("rsq", "f", "ser", "adj.rsq"),
                star.cutoffs = NULL,
                label = "tab:H2models_II",
                # star.char = c("", "", ""),
                add.lines = list(
                  c("Demographic Controls",  "No", "No", "Yes", "Yes")),
                ci=FALSE,
                omit.table.layout = NULL,
                notes.append = FALSE, 
                notes.align = "l",
                notes = c(""))

note.latex <- "\\multicolumn{5}{c} {\\parbox[t]{12cm}{ \\textit{Notes:} Demographic controls include: sex, age, race,
ethnicity, education, ideology, and region.
* $p < .1$, ** $p < .05$, *** $p < .01$}} \\\\"
tab1[grepl("Note",tab1)] <- note.latex
cat(tab1, sep = "\n")
write(tab1, '../Appendix/Tables/S2/B8.tex')

#Create table for H3
tab1<-stargazer(H3a, H3atarget, H3_conta, H3a_int,
                out= "../Appendix/Tables/S2/B9.tex",  style="qje",
                title="Treatment effects on opposing torture (H3)",
                keep = c("treat_shame", "treat_target", "treat_shame*treat_target"),
                dep.var.labels = c("Oppose Torture"),
                covariate.labels = c("Shaming", "Ally", "Shaming*Ally"),
                omit.stat = c("rsq", "f", "ser", "adj.rsq"),
                star.cutoffs = NULL,
                label = "tab:H3models_II",
                add.lines = list(
                  c("Demographics",  "No", "No", "Yes", "Yes")),
                ci=FALSE,
                omit.table.layout = NULL,
                notes.append = FALSE, 
                notes.align = "l",
                notes = c(""))

note.latex <- "\\multicolumn{5}{c} {\\parbox[t]{12cm}{ \\textit{Notes:} Demographic controls include: sex, age, race,
ethnicity, education, ideology, and region.
* $p < .1$, ** $p < .05$, *** $p < .01$}} \\\\"
tab1[grepl("Note",tab1)] <- note.latex
cat(tab1, sep = "\n")
write(tab1, '../Appendix/Tables/S2/B9.tex')

#Create table for internatinal law outcome
tab1<-stargazer(H3b, H3btarget, H3_contb, H3b_int,
                out= "../Appendix/Tables/S2/B10.tex",  style="qje",
                title="Treatment effects on international law (pre registered)",
                keep = c("treat_shame", "treat_target", "treat_shame*treat_target"),
                dep.var.labels = c("International Legal Obligation"),
                covariate.labels = c("Shaming", "Ally", "Shaming*Ally"),
                omit.stat = c("rsq", "f", "ser", "adj.rsq"),
                star.cutoffs = NULL,
                label = "tab:H3models_IIb",
                add.lines = list(
                  c("Demographics",  "No", "No", "Yes", "Yes")),
                ci=FALSE,
                omit.table.layout = NULL,
                notes.append = FALSE, 
                notes.align = "l",
                notes = c(""))

note.latex <- "\\multicolumn{5}{c} {\\parbox[t]{12cm}{ \\textit{Notes:} Demographic controls include: sex, age, race,
ethnicity, education, ideology, and region.
* $p < .1$, ** $p < .05$, *** $p < .01$}} \\\\"
tab1[grepl("Note",tab1)] <- note.latex
cat(tab1, sep = "\n")
write(tab1, '../Appendix/Tables/S2/B10.tex')


# ----------------------------------------------------------------------
# Supplementary analyses of moral licensing
# ----------------------------------------------------------------------

morals <-
  data %>% 
  filter(treat_shame==1)%>%
  group_by(treat_target) %>%
  do(tidy(lm_robust(morality ~ 1, data = .)))%>%
  mutate(.,
         Target = case_when(treat_target==0~"Shaming \n Non-ally",
                            treat_target==1~"Shaming \n Ally"))%>%
  drop_na()

ggplot(morals, aes(Target, estimate)) +
  geom_point(data = morals, size = 1) +
  geom_errorbar(data = morals,
                aes(ymin = conf.low, ymax = conf.high),
                width = 0) +
  geom_text(aes(label= round(estimate, digits = 1)), hjust=-.6, size = 3, color = "firebrick3",
            family = "Times") +
  geom_text(aes(label= paste("(",round(std.error, digits = 1),")")), hjust=-.2, vjust =2 ,  size = 3, color = "firebrick3",
            family = "Times") +
  theme_bw()+
  theme(text = element_text(size = 12, family = "Times"),
        strip.text.x = element_text(size = 12, family = "Times"),
        strip.background = element_blank(),
        axis.title.x = element_blank(),
        axis.text.x = element_text(size = 12),
        plot.caption = element_text(size = 12, family = "Times",hjust = -.02),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()) +
  ylab("U.S. Morality")

ggsave("../Appendix/Figures/S2/B8.pdf", width = 6, height = 4)


#create Table B11
a1<-lm(morality~treat_target,data)
a2<-lm(morality~treat_target + as.factor(gender) + as.factor(age2) 
       + as.factor(ethnicity) + as.factor(education) + as.factor(hispanic) +
         as.factor(political_party) + as.factor(region), data = data)

tab1<-stargazer(a1, a2,
                out= "../Appendix/Tables/S2/B11.tex",  style="qje",
                title="The Effect of Ally-shaming on Perceptions of U.S. Morality",
                keep = c("treat_target"),
                dep.var.labels = "Government Approval",
                covariate.labels = c("Ally"),
                omit.stat = c("rsq", "f", "ser", "adj.rsq"),
                star.cutoffs = NULL,
                label = "tab:morals_target",
                add.lines = list(
                  c("Demographic Controls",  "No", "Yes")),
                ci=FALSE,
                omit.table.layout = NULL,
                notes.append = FALSE, 
                notes.align = "l",
                notes = c(""))

note.latex <- "\\multicolumn{3}{c} {\\parbox[t]{12cm}{ \\textit{Notes:} Demographic controls include: sex, age, race,
ethnicity, education, ideology, and region.
* $p < .1$, ** $p < .05$, *** $p < .01$}} \\\\"
tab1[grepl("Note",tab1)] <- note.latex
cat(tab1, sep = "\n")
write(tab1, '../Appendix/Tables/S2/B11.tex')

#Mediation (an exploration)
#create figure B10

morality_med <- lm(morality~treat_shame, 
                   data=data)
data$age<-as.numeric(data$age2)
data$torture<-as.numeric(data$torture)
morality_out <- lm(torture~treat_shame + treat_target+ morality + gender + political_party +
                     education + region, data = data)

mediation_morality <- mediate(morality_med, morality_out, treat = "treat_shame", 
                              mediator = "morality", dropobs=TRUE, 
                              boot=TRUE, conf.level=.90)
summary(mediation_morality)

#save pdf
pdf(file="../Appendix/Figures/S2/B10.pdf")
plot(mediation_morality, labels=c("ACME\n(Morality)",
                                  "Direct \nEffect", "Total \nEffect"), xlim = c(-1,1))
dev.off()

##Sensitivity analysis
sens.out <- medsens(mediation_morality, rho.by = 0.1, effect.type = "indirect", sims = 100) 
par(mar=c(5,6.5,3.2,01))

plot(sens.out, sign.prod="positive", sens.par = "R2", r.type = c("residual"), 
     xlab="Proportion of Total Variance in M Explained by Confounder", 
     ylab="Proportion of Total Variance in Y Explained by Confounder", main="(a)")

plot(sens.out, sign.prod="negative", sens.par = "R2", r.type = c("residual"), 
     xlab="Proportion of Total Variance in M Explained by Confounder", 
     ylab="Proportion of Total Variance in Y Explained by Confounder", main="(b)")


# ----------------------------------------------------------------------
# Manipulation Checks
# ----------------------------------------------------------------------

#Did the U.S. government criticize the other country?
data$manipulation_1 = recode(data$attention_2, "1=3; 2=2; 3=1") #recode such that [yes(3), not specified (2), no (1)]

man<- lm(manipulation_1~treat_shame, data = data) 
summary(man)

man_controls<-lm(manipulation_1~treat_shame+ as.factor(gender) + as.factor(age2) 
                          + as.factor(ethnicity) + as.factor(education) + as.factor(hispanic) +
                            as.factor(political_party) + as.factor(region),data=data)

summary(man_controls)

#Was the other country a US ally?
data$manipulation_2 = recode(data$attention_3, "1=3; 2=2; 3=1") #recode such that [yes(3), not specified (2), no (1)]

man_target<- lm(manipulation_2~treat_target, data = data) 
summary(man_target)

man2_controls<-lm(manipulation_2~treat_target+ as.factor(gender) + as.factor(age2) 
                          + as.factor(ethnicity) + as.factor(education) + as.factor(hispanic) +
                            as.factor(political_party) + as.factor(region),data=data)

summary(man2_controls)

tab1<-stargazer(man, man_controls,
                man_target, man2_controls,
                out= "../Appendix/Tables/S2/B12.tex",  style="qje",
                title="Manipulation Check",
                keep = c("treat_shame", "treat_target"),
                dep.var.labels = c("US criticized the other country",
                                   "Other country is a US ally"),
                covariate.labels = c("Shaming treatment", "Ally treatment"),
                omit.stat = c("rsq", "f", "ser", "adj.rsq"),
                star.cutoffs = NULL,
                label = "tab:manipulationII",
               # star.char = c("", "", ""),
                ci=FALSE,
                add.lines = list(
                  c("Demographic Controls",  "No", "Yes","No", "Yes")),
                omit.table.layout = NULL,
                notes.append = FALSE, 
                notes.align = "l",
                notes = c(""))

note.latex <- "\\multicolumn{5}{c} {\\parbox[t]{12cm}{ \\textit{Notes:} Demographic controls include: sex, age, race,
ethnicity, education, ideology, and region.
* $p < .1$, ** $p < .05$, *** $p < .01$}} \\\\"
tab1[grepl("Note",tab1)] <- note.latex
cat(tab1, sep = "\n")                
write(tab1, '../Appendix/Tables/S2/B12.tex')

#how many answered manipulation checks correctly?
data<-data%>%
  mutate(.,
         pass_shame = case_when(attention_2=="1" & treat_shame == "1"~"1",
                          attention_2=="2" & treat_shame == "1"~"0",
                          attention_2=="3" & treat_shame == "1"~"0",
                          attention_2=="1" & treat_shame == "0"~"0",
                          attention_2=="2" & treat_shame == "0"~"0",
                          attention_2=="3" & treat_shame == "0"~"1"),
         pass_ally = case_when(attention_3=="1" & treat_target == "1"~"1",
                               attention_3=="2" & treat_target == "1"~"0",
                               attention_3=="3" & treat_target == "1"~"0",
                               attention_3=="1" & treat_target == "0"~"0",
                               attention_3=="2" & treat_target == "0"~"0",
                               attention_3=="3" & treat_target == "0"~"1"))

sum(as.numeric(data$pass_shame),na.rm=T)/sum(!is.na(data$pass_shame))
sum(as.numeric(data$pass_ally),na.rm=T)/sum(!is.na(data$pass_ally))

# ----------------------------------------------------------------------
# Attrition
# ----------------------------------------------------------------------

#data1 indicates all respondents who passed attention. 
data1 <- data1 %>%
  mutate(.,
         treat_shame = case_when(
           condition == "no_shame_ally" ~ 0,
           condition == "no_shame_no_ally" ~ 0,
           condition == "shame_ally" ~ 1,
           condition == "shame_no_ally" ~ 1),
         treat_target = case_when(
           condition == "no_shame_ally" ~ 1,
           condition == "no_shame_no_ally" ~ 0,
           condition == "shame_ally" ~ 1,
           condition == "shame_no_ally" ~ 0),
         attrite = ifelse(is.na(data1$religion),1,0)) #religion is the last question in the survey

sum(data1$attrite) #146 respondents
(sum(data1$attrite)/nrow(data1))*100 #only 4.6% of the sample

#Does my treatment cause people to attrite?
attrite_model<-lm(attrite~treat_shame,data=data1)
summary(attrite_model) #no

attrite_controls<-lm(attrite~treat_shame+ as.factor(gender) + as.factor(age2) + as.factor(ethnicity) + as.factor(hispanic) +
                       as.factor(political_party) + as.factor(region),data=data1)
summary(attrite_controls) #also when controling for demographics

#Create table
tab1<-stargazer(attrite_model, attrite_controls,
                out= "../Appendix/Tables/S2/B13.tex",  style="qje",
                title="The effect of shaming treatment on attrition",
                keep = "treat_shame",
                dep.var.labels = "Attrition",
                covariate.labels = "Shaming treatment",
                omit.stat = c("rsq", "f", "ser", "adj.rsq"),
                star.cutoffs = NULL,
                label = "tab:attriteII",
                # star.char = c("", "", ""),
                ci=FALSE,
                add.lines = list(
                  c("Demographic Controls",  "No", "Yes")),
                omit.table.layout = NULL,
                notes.append = FALSE, 
                notes.align = "l",
                notes = c(""))

note.latex <- "\\multicolumn{3}{c} {\\parbox[t]{12cm}{ \\textit{Notes:} Demographic controls include: sex, age, race,
ethnicity, education, ideology, and region.
* $p < .1$, ** $p < .05$, *** $p < .01$}} \\\\"
tab1[grepl("Note",tab1)] <- note.latex
cat(tab1, sep = "\n")                
write(tab1, '../Appendix/Tables/S2/B13.tex')

# ----------------------------------------------------------------------
# Exploring confounding
# ----------------------------------------------------------------------

data1$placebo = ifelse(!is.na(data1$placebo_1_TEXT),1,0)

counfoudning<-lm(placebo~treat_shame, data=data1)
summary(counfoudning)

counfoudning_controls<-lm(placebo~treat_shame+ as.factor(gender) + as.factor(age2) 
                     + as.factor(ethnicity) + as.factor(education) + as.factor(hispanic) +
                       as.factor(political_party) + as.factor(region),data=data1)

summary(counfoudning_controls) 

#Create table
tab1<-stargazer(counfoudning, counfoudning_controls,
                out= "../Appendix/Tables/S2/B14.tex",  style="qje",
                title="The effect of shaming treatment on thinking of a specific country",
                keep = "treat_shame",
                dep.var.labels = "Did you think of a specific country?",
                covariate.labels = "Shaming treatment",
                omit.stat = c("rsq", "f", "ser", "adj.rsq"),
                star.cutoffs = NULL,
                label = "tab:confoundingII",
                #star.char = c("", "", ""),
                ci=FALSE,
                add.lines = list(
                  c("Demographic Controls",  "No", "Yes")),
                omit.table.layout = NULL,
                notes.append = FALSE, 
                notes.align = "l",
                notes = c(""))
note.latex <- "\\multicolumn{3}{c} {\\parbox[t]{12cm}{ \\textit{Notes:} Demographic controls include: sex, age, race,
ethnicity, education, ideology, and region.
* $p < .1$, ** $p < .05$, *** $p < .01$}} \\\\"
tab1[grepl("Note",tab1)] <- note.latex
cat(tab1, sep = "\n")                
write(tab1, '../Appendix/Tables/S2/B14.tex')
